function [wdraw,what_star,wplus]=simuldis(y,Z,T,R,Qmat);
% ==================================================================
% This function takes as inputs the matrices of a state space model with no
% error in the observation equation and the Kalman Gain, and F^(-1) matrices from 
% a forward run of the Kalman Filter in order to simulate a set of draws
% from the state disturbances using a disturbance smoother. 
%
% The models is 
% y(t) = Z*s(t) 
% s(t) = T*s(t) + R*n(t)  V( n(t) ) = Q(t)  
%
% Where y is NZ x 1, x is NY x 1 and n is NX x 1 
% 
% Inputs 
% ======
% y     nobs x NZ stacked vector of observations, nobs is the total number of observations  
% 
% All matrices but Q(t) are assumed to be time invariant to simplify the computation 
% The correct size of input Qmat then depends on whether Q is time varyong or  not. 
% not time varying    Qmat is NX x NX 
%      time varying   Qmat is NX x NX x T where
% 
% Alejandro Justiniano and Giorgio Primiceri (c) April 2005 
% =============================================================
[nobs,nz]=size(y); 
[ny,nx] = size(R); 

ch=size( Qmat ); 
if ch(2) ~= nx; error('Column dimension of Qmat must be nx'); end; 
if ch(1) == nx      
    flag_tv = 0 ;
elseif ch(1) == nobs   
    flag_tv = 1; 
else 
    error('Row dimension of Qmat must be nx if not TV or nobs if TV') 
end 
clear ch; 

% =========================================
% Simulate the innovations and the data
if flag_tv == 0 
        wplus=randn(nobs,nx)*Qmat; 
else 
        wplus=randn(nobs,nx).*Qmat; 
end 
clear qtemp temp; 


% This is y+ in DK. Note that the first state is drawn from a 
% Normal ( 0, eye(2) ) 

ny=size(T,1); 
[yplus,stplus]=simssmodel(wplus,T,R,Z,zeros(ny,1),eye(ny));

% According to the paper, if a single draw is needed, 
% ystar = y - yplus generated above 
ystar=y-yplus; 

% =================================================================
% Run the forwrad recursion of the filter and obtain the matrices
[junk,vstar,finmat,kmat]=kfilter_mod( ystar , Z, T, R, sqrt(2)*Qmat ); 

% ==================================================================
% Compute the mean of the disturbances using the backward smoother 

what_star = zeros( nx,nobs); 

Ztr=Z'; %clear Z; 
Rtr=R'; %clear R; 
Ttr=T'; %clear T; 

rstar=zeros(ny,1); 

ii=nobs;

switch flag_tv 
    case 0;        
        Q=Qmat'*Qmat; 
        for ii=nobs:-1:1;  
            Ltr=Ttr-Ztr*kmat(:,:,ii)';
            [rstar,what_star(:,ii)] = smoothdis( rstar, Q, Rtr, Ztr, finmat(:,:,ii), Ltr , vstar(:,ii) ); 
        end 
    case 1; 
        for ii=nobs:-1:1;  
            Ltr=Ttr-Ztr*kmat(:,:,ii)';    
            J=2*diag( Qmat(ii,:).*Qmat(ii,:) );
            [rstar,what_star(:,ii)] = smoothdis( rstar, J, Rtr, Ztr, finmat(:,:,ii), Ltr , vstar(:,ii) ); 
        end 
end
wdraw = wplus + what_star' ; 